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We study the charmonium spectrum including higher spin and exotic states. We use the 
Sheikholeslami- Wilson (clover) action for Nf — 2 sea quarks as well as for the charm valence 
quark. In order to access excited states we apply a variational method with a basis of highly 
optimized operators. 
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1. Introduction 

The discovery of about ten new charmonium resonances within the past five years prompted 
an increased experimental and theoretical interest in the phenomenology of these states. This trend 
is bound to continue, in particular since the PANDA experiment at the new antiproton facility 
FAIR will produce huge new data samples whose interpretation will require theoretical input. One 
way both to reproduce and to predict the experimental spectrum of particle resonances from first 
principles QCD are Lattice simulations. 

Properties of many recently observed states are at variance with nonrelativistic quark model 
predictions, the most striking example being the X(3872). It is an important task to reveal the inner 
structure of such states, i.e. to clarify to what extent charmonium resonances can be interpreted 
as quark model cc states or whether some of these contain significant quark-gluon hybrid or four 
quark contributions (tetraquark or molecule). 

We expect the lightest hybrid charmonium states to be heavier than 4.3 GeV, the experimental 
DD threshold lies above 3.7 GeV and the vector charmonium J/y ground state much lower, at 
about 3.1 GeV. This means that, with the exception of D waves and higher angular momentum 
states, the higher Fock components will only start to show up prominently in radial excitations. As 
a first step, we need to be able to reliably compute these on the lattice. 

We describe the methods used to obtain these states, our simulation set-up and present first 
results on the charmonium spectrum and "wavef unctions". 

2. Variational method 

We start from a cross correlator matrix 

Qjit) = {0i{t)0]m, (2.1) 

with a basis of operators Oi, i = I, . . . ,N, destroying a colour singlet state within the desired lattice 
Oh^C representation from which we wish to deduce continuum quantum numbers. We do not 
include charmed sea quarks and hence the numbers of charm and anti-charm quarks are separately 
conserved. At present we restrict ourselves to the cc sector (including hybrids). At a later stage 
four quark operators will be incorporated. 

The correlator matrix can be spectrally decomposed, 

Cij{t)='£vy;e-^"', (2.2) 

n 

where v" is the nth state within the subsector of the Hilbert space spanned by C{t), and £"„ is 
the corresponding energy eigenvalue. Since Ctj is a real symmetric matrix^ the v" are mutually 
orthogonal. 

To obtain initial guesses of the eigenvalues A" and -vectors a = I,. . . ,N, we solve the 
generalized eigenvalue problem ^ ^, 

C-'/^fa)C(O C-'/^(foX(0=A«(0<(0, (2.3) 

'This only holds in the limit of infinite statistics. We symmetrize C{t) by hand, after checking that violations are 
consistent with zero, within the statistical errors. 
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varying Iq and t > Iq. Due to this symmetrized construction the A/^-component eigenvectors 
mutually orthogonal for all choices of to and t. Note that the eigenvalues of the system C(f)(^" = 
A"C(fo)^" with = C'/^(fo)0" are the same as in Eq. ( |23| ). Moreover, the non-orthogonal ^s 
will approach the orthogonal i/As at large Iq. 

If we choose to too large, the rank of C{to) will not be maximal anymore as (within statistical 
errors) excited states will die out in Euclidean time. For to chosen too small, C{t) will receive 
contributions from more than the A'^ lowest lying states, resulting in unstable eigenvectors and 
effective masses, 

<«<'>=°"'»(x^)' 

where we have suppressed the subscript to- For each channel we employ a three dimensional basis 
of operators which we call local, narrow and wide. 

We apply iterative Gaussian smearing to the fermion fields (j), 

0i"^''=c(^0i"' + '^.f^.vA¥^^^ (2-5) 

with K = 0.3 and a normalisation c to avoid numerical overflow. We smear quark and antiquark with 
the same number of steps ng which is equivalent to applying 2ng smearing steps to one propagator 
only. 

The parallel transporters Uxj = uj^^p above are p-APE smeared: 



where a = 2.5 was chosen to maximize the spatial plaquette constructed from the smeared links 
(see the left hand side of Figure [T]). Psu{3) denotes a projection operator, back into the gauge group. 



By using smeared transporters within Eq. (2.5) we achieve a more continuum-like spatial distri- 
bution of the smearing wavefunction (see Section |5] below) and better overlaps with the physical 
states. 

Effective masses for symmetric 2-point functions with a local source and different sink smear- 
ings are shown in Figure |l]. We selected the smearing applied to the trial wavefunctions within 
each channel, such that one effective mass approached the asymptotic state from above, one from 
below and ones sat exactly on spot. This procedure ensures that the span of our variational basis 
has overlap not only with the ground state but also with the lowest radial excitations. 



3. Simulation details 

While chiral symmetry plays a minor role for charmonia, the charm quark mass is not 
heavy enough to allow for a non-relativistic treatment. We use the clover Wilson action both for 
valence and sea quarks which will give us a well-defined continuum limit. However, mc-a <^ 1, such 
that ultimately an a — > extrapolation will be important. We work on Nf = 2 dynamical lattices 
generated by the QCDSF collaboration. Details of these lattices can be found in Ref. Here 
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Figure 1: Dependence of the average value of the spatial plaquette on the number of p-APE smearing steps 
(left) and effective masses for local-smeared correlators with different smearing functions (right). 



we present results from a single 16^ x 32 lattice at j8 = 5.20 and K = 0.13420 corresponding to an 
inverse lattice spacing of a"' w 1.73 GeV and a pion mass rriji ~ 1 GeV. The lattice spacing was 
determined from the value ro ~ 0.46 fm such that the nucleon reaches its experimental mass when 
extrapolated to physical mj^. This leaves us with the charm quark mass as the only free parameter 
which we set by tuning mj^ = \m^^ + \mjj^i to the experimental value. 

The operators we use are based on Ref. [Q], however derivatives were symmetrized to allow for 
charge conjugation eigenstates also at finite momenta. The quark bilinears about which we report 
here are displayed in Table |l], together with their irreducible lattice representations and the lowest 
spin continuum state they couple to. 
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Table 1: Interpolating fields in use. 



4. Spectrum 

In Figure ^ we show effective masses for the analyzed lattice. In all channels we see nice 
plateaus for the ground and the first excited state. Apart from the exotic 1^+ channel we also 
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obtain reasonable signals for the second excited state which however we digest with caution: to 
gain more confidence in these we will move to a larger operator basis and increase statistics. The 
lines indicate the fit ranges and errors. For the 1^+ it turned out particularly hard to separate the 
ground state from the first excitation because these are very close in mass. Hybrid potentials are 
rather flat and yield dense spectra within potential models. Hence this maybe taken as a hint at a 
hybrid content. However, the creation operator with best overlap with this state does not contain 
an explicit chromomagnetic field dependence: in a relativistic theory, for any allowed J all PC 
quantum numbers including exotic ones can be obtained from quark bilinears, even in the free field 
case. More study of this question is required. 

The computed spectrum is plotted in Fig. |3|, together with the experimental values. It is im- 
portant to note that we were not too careful when setting the charm quark mass parameter and 
underestimate m-^ by about 15 MeV. However, this is still well within the accuracy of the lat- 
tice spacing determination and moreover will cancel from level splittings. So the whole spectrum 
should be shifted higher a bit. 

We obtain a IS hyperfine splitting of Amis = 73(2) MeV, below the experimental 117 MeV 
and consistent with the unphysically high pion mass, wrong number of sea quarks and lattice arte- 
facts. Also disconnected contributions can affect this quantity. For the 2S hyperfine splitting we 
obtain Am2s = 47(6) MeV, in agreement with experiment [49(5) MeV]. Our failure to consistently 
underestimate this value as well might be explainable by DD threshold effects which we neglect 
due to our heavy sea quarks (and hence D mesons). Alternatively, the rjc mass might receive small 
contributions from the L'a(I) anomaly or from rjc — rj' mixing, also effects that we neglect, whose 
full treatment requires disconnected quark line diagrams and charmed sea quarks. 
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Figure 2: Effective masses from the three dimensional operator basis. Fit ranges and errors are indicated by 



horizontal lines. The to values refer to the respective normaUzation time slices (see Eq. (2.3)) 
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DD 



Figure 3: Predicted spectrum, together with the experimental values. The DD threshold is the experimental 
one. 

5. Wavefunctions 



The variational method not only helps to compute the spectrum but also provides access to 
couplings. Efforts to extract such couplings in the heavy quark regime have already been made 
for example in Refs. and [^]. There exists no real shortcut to the computation of three-point 
functions for this purpose but nonetheless it can be instructive to analyze the smearing functions 
that we use in some detail. 



For sufficiently large to we can identify the components of the eigenvectors of Eq. (2.3) as 
the couplings of our interpolating functions with the physical state of mass m^. Here we start from a 
four-dimensional basis of trial functions with 0, 5, 10 and 40 (times two) Gauss smearing iterations. 
We can apply these functions to 5 sources (in space and colour) to obtain the spatial distribution 
of the corresponding trial wavefunctions <I>;(x) € SU{3). By folding these with our eigenvectors 
we can attempt to construct the "wavefunctions" of the physical states: *P"(x) = i/a"<I>j(x). 
Needless to say that it is not possible to exactly create the physical eigenstates with such a small 
number of trial functions. Moreover, in Euclidean time we cannot obtain the phase information and 
hence only the (gauge invariant) probability densities are meaningful quantities. Thus in general 
Xr^p^txj//? ^ 5a/3 Unfortunately, the used lattice is too coarse to resolve the node structure of 
However, observing that for our APE smeared fields and in the free case are very similar, 
we plot the wavefunctions for the free case, where the nodes are clearly visible due to the sign 
change. Fig. ^ shows these wavefunctions for the lowest three pseudoscalar states. We neglect the 
statistical errors. 

In spite of the small basis the node structure is consistent with the IS, 2S and 3S assignment, 
with no visible pollution from higher Fock states or D waves. For the 15 we obtain an rms width 
of ca. 0.39 fm. This compares reasonably well with the infinite volume continuum potential model 
expectation of about 0.4 fm 
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Figure 4: The IS, 2S and 3S pseudoscalar "wavefunctions". 



6. Outlook 

The next step is to include disconnected contributions to the charmonium two-point functions 
which may play an important role for quantities like the IS hyperfine splitting. Work on improved 
all-to-all propagator calculations is already in progress. We also plan to extend our operator basis 
by including four-quark states, once we have moved to light sea quark masses. 
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